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Abstract. The radiative transfer equation (RTE) for po- 
larized light accepts a convenient exponential solution 
when the absorption matrix commutes with its integral. 
We characterize some of the matrix depth variations which 
are compatible with the commutation condition. Eventu- 
ally the vector solution may be diagonalized and one may 
obtain four independent scalar solutions with four optical 
depths, complex in general. When the commutation con- 
dition is not satisfied, one must resort to a determination 
of an appropriate evolution operator, which is shown to be 
well determined mathematically, but whose explicit form 
is, in general, not easy to apply in a numerical code. How- 
ever, we propose here an approach to solve a general case 
not satisfying the commutation condition. 

Key words: Sun:magnetic fields - Techniques: polarimet- 
ric - Radiative transfer - Polarization - Methods: numer- 
ical 



1. Introduction 

The use of spectropolarimetry and Zeeman effect to mea- 
sure magnetic fields in the sun and stars, requires a trans- 
fer theory for polarized light, in the presence of a magnetic 
field. Such transfer equation was first written in the pio- 
neering paper by Unno ( 1956 ), where the effect of anoma- 
lous dispersion was not considered yet. Rachkowsky ( |1962| ; 



1967) was the first to include it, and obtained a general 



transfer equation for polarized light. 



— / = -K/+ J, 
dz 



(1) 



in terms of the Stokes parameters I, Q, U and V presented 
as the components of vector /. This equation resembles 
the transfer equation for the scalar case when polarization 
is ignored. The scalar intensity is replaced by the Stokes 
vector I. The emision term is now a four component vector 
J, and the scalar absorption coefficient becomes a 4 x 4 
matrix K which describes absorption (including anomalous 
dispersion) in the presence of Zeeman effect. The variable 
z parameterizes the light path. The transformation from 



a single equation to a system of four equations (abridged 
in a vectorial form) changes drastically the nature of the 
problem and no general explicit analytical solution has 
been proposed so far. 

The first particular solution was obtained by Unno, 
who applied the equation to the case of a Milne-Eddington 
atmosphere. Rachkowsky (1967), after including anoma- 
lous dispersion, solved the RTE for the same homogeneous 
atmosphere. 

Several procedures were successful in solving the 
equation numerically (e.g. Beckers & Schroter 1969, 
Wittman 1974, Rees, Murphy & Durrant 1989, Landi 
DeglTnnocenti 1976) even for the general case, by subdi- 
viding the atmosphere into numerous layers. The lasts are 
chosen optically thin in the upper atmosphere and eventu- 
ally an Unno-Rachkowsky solution is taken for the deep- 
est layer, down to optical depth infinity. Difficulties were 
encountered as well, depending on the particular method 
chosen to solve the radiative transfer in each individual 
sublayer. An universal technique, like the Runge-Kutta 
method, may not be the best approach for a particular 
case, mainly because of different scales of variation. The 
mathematical justifications are often not rigourous and 
numerical tests are always necessary. 

The advantage of an eventual analytical solution is 
obvious. But, up to now, they have been always re- 
stricted to homogeneous atmospheric models where only 
the LTE source function was depth dependent. The con- 
stant K matrix has been handled with differe nt m athe- 
matical techniques : for instance Kjeldseth Moe ( |l968| ) and 
Ste nflo(1971 ; 1994) used K-diagonalization, van Ballegooi- 
jen( 1985 ) preferred Jones calculus. 

Now, for the solar case, variations with depth of both 
the thermodynamical parameters describing the atmos- 
phere and the magnetic field, are not negligible. For in- 
stance, Ruiz Cobo & del Toro Iniesta( i992| ; 1994), and del 
Toro Iniesta & Ruiz Cobo (199£ ; 1996D using numerical in- 
version of the observed Stokes profiles have confirmed the 
need for inhomogeneous models (see Collados et al. 1994, 
Westendorp Plaza et al. 1997a, 1997b, 1997c, and see also 
del Toro Iniesta & Ruiz Cobo 1996a for a review). All 
these works raise the interest in analytical methods deal- 
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ing with non-homogeneous atmospheric models, i.e. with 
non-constant K matrices. 

At first sight, the RTE for polarized light does not 
seem more complicated to solve than its scalar equivalent 
(RTE for pure light intensity, with no polarisation): 



dz 



I = -kI + J, 



where as usual, k is the absorption coefficient, / and J 
are the scalar intensity and emission function respectively 
and z the geometrical path. This equation has an explicit 
formal solution, for the layer zq < t < z : 



I{z) = 



— \ ndt' 
3 Jt 



J{t)dt + e 



— f Kdt' 
•1-0 I 



(2) 



A direct extrapolation of this scalar solution would yield: 
I{z)= re-r^''*'j(t)dt + e"^^o^''*/(zo). (3) 

Note that, in this expression, all the mathematical opera- 
tions involving matrices are well defined. For instance, the 
exponential of a square matrix is defined as 



(4) 



Now, just one difference arises between algebras using ma- 
trices and scalars: while two scalars always commute, that 
is not true in general for two matrices: 

KL LK. 

This difference becomes important when treating for ex- 
ample the derivative of a power of a matrix: it is not true 
in general that 



dz 



L" = -nKL" 



where we have used (and so we will do hereafter) that 
d 



dz 



L = -K. 



As a further consequence of the non-commutativity of ma- 
trices, we have that 



dz 



in complete contradiction with the scalar case. 

In short, Eq.(3) is a solution of the polarized RTE, 
Eq.(l), when the commutation condition. 



[K,L] = KL - LK = 0, 



(5) 



holds. Under this assumption, the previous expressions re- 
cover an usual scalar appearance, and therefore we are per- 
mitted to write a scalar-like formal solution. It is interest- 
ing to note that this condition docs not imply a constant 
absorption matrix. In the following sections wc will show 



how to incorporate variations of K with optical depth. In- 



deed, Landi DeglTnnoccnti & Landi deglTnnocenti (1981 



1985) have already shown how to handle matrices of the 
K = K'f{z) 



form 



with K' being a constant matrix. Although not explicitly 
said in these papers, it is obvious that here K satisfies 
condition (||). 

For the more general case when the commutation con- 
dition (||) does not hold, Landi deglTnnocenti & Landi 
deglTnnocenti (198£) have derived a formal solution for 
the RTE: 



/(z)= / 0{z,z')J{z')dz' + 0{z,zo)I{zo), 



(6) 



where 0{z, z') is the evolution operator, a new 4x4 matrix 
which obeys the homogeneous equation 



— O(z,zo) = -K(z)0(z,zo), 
dz 

with initial condition 



(7) 



0(zo,2o) = 1, 

where 11 represents the 4x4 identity matrix. Note that 
when solution (||) applies, the evolution operator takes an 
explicit form, namely: 



0(z, z') = exp(L(z, z')) = exp 



Kdt 



(8) 



A general method to solve equations of the form of Eq. 
( 7\j for linear operators has already been given by Magnus 
(1954 ) . In this remarkable paper an exponential expression 
is proposed: 

0(z, z') = exp {n{z — z')) , 



were the exponent Q{z — z') is given by an infinite series. 

Equation (|^) turns out to be Magnus expression when 
only the first term in the infinite series is kept (indeed the 
only non-zero one when condition (^ holds). 

We now discuss three existing options to solve Eq.(n]): 



Constant matrix assumption. Condition (5) is inmedi- 
ately satisfied and analytical solutions were found 
( |Unno 19561 ; [Rachkowsky 1967| ). 

Multi-layer techniques. The atmosphere is considered 
as made up by numerous successive layers. A crude 
assumption on the radiative transfer in each optically 
thin layer is then advanced, and leads to a procedure 
of numerical integration expected by intuition to con- 
verge to the exact solution. A formal proof of con- 
vergence was not given, but the numerical tests were 
indee d sati sfactory. See for instance Rees ( 1987 ), Rees 
et a?.( |l989[) , Ruiz Cobo & del To ro Iniesta7|l99|) , del 
Toro Inicsta & Ruiz Cobo (|l996| ). 
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3. Magnus solution. By applying linear algebra one can 
treat the general case, with non commuting matrices 
([Magnus 1954|). 



The constant matrix technique, method (1), is not possible 
when one wants to abandon the homogeneous magnetic 
field and atmosphere assumptions. Next, poor economy is 
the main drawback of method (2). There is some doubt 
whether one can determine a priori the number of layers 
necessary for a desired precision. Last but not least, a more 
analytical insight than what a pure numerical method can 
give is always desired as well. Moreover, our ultimate pur- 
pose is magnetometry of the sun or stars. We want to 
go beyond the first method, the most used, at present, 
but limited to a constant absorption matrix and therefore 
also constant field. Still we must admit that actual obser- 
vations will allow us to determine only little more than 
a homogeneous atmosphere model, say, at most the mag- 
netic fields at two or three levels in the atmosphere. It 
is therefore not "economic" to calculate more than a few 
layers in the atmosphere. 

It is striking that Magnus' solution was published two 
years before the memorial paper by Unno(1956), the first 
paper on RTE for polarized light, and as yet it has never 
been mentioned in the astrophysical literature. It is there- 
fore given in Appendix A. The solution given by Magnus 
is mathematically exact, but it requires the use of Lie al- 
gebra, is not economic and can hardly be used in practical 
computation. It is mentioned here because it confirms the 
approach of the present paper and complete it. 

Our general strategy is, first, to "satisfy" condition (5) 
as far as possible by extracting from the absorption matrix 
everything that commutes with its integral and therefore 
can easily be integrated according to Eq.(3), as explained 
in section 2. In section 3, we diagonalise the commutative 
part of the matrix to allow an efficient integration. Then, 
in section 4, we treat the residual matrix by an appropriate 
approximation and thus obtain a semi-analytical solution 
for an optically finite layer with arbitrary depth varia- 
tions. Eventually we can then borrow the techniques from 
the multi-layer approach and apply our semi-analytical so- 
lution to a few layer model to improve the computation. 

A few words on the mathematical space where we are 
working and where the RTE is to be solved, are in order. 
Magnetometry concerns the 3D real physical space, where 
the magnetic field can be represented as a 3D "vector" 
and all physical parameters of the atmosphere determine 
the coefficients that enter the radiative transfer equation. 
The last one is much better calculated in another space. 
Indeed, we have already entered another 4D geometry: the 
Minkowski space, where the Stokes' 4-vectors are best de- 
scribed. In this geometry, the norm of a vector I is given 
by [P — ~ ~ V^)^ ■ It has particular symmetries 
and is governed by linear algebra. The elementary oper- 
ations, like absorption and retardation, are presented by 
matrices for which commutation relations are of particu- 



lar importance. When condition (5) holds, an exponential 
solution, scalar like, to a linear equation can easily be de- 
rived. Otherwise, we have to turn to Magnus' exponential 
solution. 

The main difficulties originate from the fact that only 
few variables are explicitly common to the "two spaces" . 
Typically scalar variables like z,kc,Ki (see section 2. for 
their definitions) will appear in both spaces in similar 
ways. However, rotations of the Stokes reference system 
will not. Exception is the azimuth rotation. The angle of 
rotation of the azimuth of the magnetic field in the "real 
3D space" corresponds to a rotation in the Minkowski 
space, but with a double amount. Naturally, when a con- 
stant atmosphere is selected in the real space, the cor- 
responding matrix in the Minkowski space will be con- 
stant as well. On the other hand some rotation in the 
Minkowski space may be much easier to handle. For in- 
stance, one may find convenient to use generalized Stokes 
vectors expressed in terms of elliptic states of polariza- 
tion. Transformations from one set of Stokes reprensenta- 
tion to another are expressed simply as rotations in the 
Minkowski space. Except in some limiting cases it is not 
possible to translate these angles in terms of angles in the 
physical space. At the same time, the highly non linear 
relations between magnetic field and the entries of the 
absorption matrix cannot in general be simplified. Thus, 
while the RTE can be solved for a given depth variation 
of the absorption matrix, we cannot, in general, recover 
analytically the corresponding variation of the magnetic 
field. We anticipate that numerical methods can overcome 
this difficulty and profite from the analytical solution in 
the Minkowski space to treat the depth variations of the 
magnetic fields and improve both the economy and the 
precision of the calculations. These considerations apply 
as well to all other atmospheric conditions, like tempera- 
ture, pressure, velocity etc. 

In some particular cases, the relations between vari- 
ables in the Minkowski and real spaces may become sim- 
plified. For instance, in absence of absorption of linear 
polarization, whether in the pure longitudinal magnetic 
field, or alternatively for particular Zeeman patterns, free 
of linear polarization. Also for the case when all Zeeman 
components are separated, simple relations hold as will be 
discussed in the corresponding sections. 

2. Transformation of matrix K 

We rewrite the transfer equation as 

d 



dz 



Iq = — Kq/q + Jo 



And, Kq being invertible, we can define the source function 
vector, either LTE or not, 

So = Kq ^ Jo, 
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so that the transfer equation reads 



so that we can write the transformed transfer equation as 



— /o = -Ko(Jo-5o). 

dz 

Matrix Kq can be decomposed as follows 
/ 



(9) 



Ko = Ki{z) 





h cos 2(f) 
b sin 20 
c 



b cos 2<j) 


-7o 
(3 sin 20 

+{gKi{z) - 



b sin 20 

7o 


-/3 cos 20 
«;,(z))]l 




where 



with 



Ki = K((z)Ni + {gKi{z) + Kc{z))l, 



Ni 



(14) 



b 





c 





7 





-7 





P 





-/3 






Where ^((-z) and Kc{z) are the usual scalar absorption 
coefficients: the selective (at line center) and the contin- 
uum one, respectively; is the azimuth angle of the mag- 
netic field, relative to a fixed reference system, and 1 is 
the 4x4 identity matrix. 

This is the general symmetry of Kq; the meaning of 
parameters Oj5, c, /3,and 70 can be found by comparing 
expression (10) with the correspondi ng ones in Landi 
DegTInnocenti & La ndi D eglTnnocenti( 1981 ; 1985 ), Rees 
( |1987| ) or Kawakami( pJ83|) . 

We can simplify this matrix by rotating it an angle 20 
in the plane Q-U. That is, we introduce a rotation matrix 



Ri 



and its inverse Rj^ and we apply them to Kq to obtain 






K'l = RiKoRr' = Ki{z) 




(11) 



Applying this transformation to Eq. (^ we obtain: 

R^/o 

az 



where 



(RiKoRr') Ri/o + Ri Jo = -K;/i + Ji (12) 

Ji = RiJo, 
Ji = Ri«/o- 

The left hand side of the transformed transfer equation 
(O) is equal to 



(Ri/o) 



dz 

where we note that 

(:^Ri)/o = (:^Ri)Rr'RiJo 

dz dz 



dz 



(;fRiUc 

dz 




(13) 



where we have introduced 



7 = 70 



2^-. 

dz Ki 



(15) 



The meaning of the new Stokes reference system is as fol- 
lows: after the rotation, the new generalized Stokes pa- 
rameters Qi and Ui, projections of vector / on axis Qi 
and Ui in the new reference system, still correspond to 
linear polarization, but Zeeman linear absorption affects 
Qi only (absorption along the Qi axis). Faraday rotation 
may still affect Ui, but with zero absorption. The parame- 
ters Ii and Vi are unchanged (the correspoding axis I and 
V are not affected by the rotation). In the real space, 
the meaning of this rotation is that the reference for the 
usual definition of the Stokes parameters is taken parallel 
to the magnetic field for Q. These new axes rotate with 
the field. 

A second simplification is obtained by the use of a new 
rotation, given by 



(16) 




and its inverse R2 , where 
6 



cos a 



sm a 



(17) 



By applying it to the matrix Ki we obtain 



K2 = R2K1R, 



Kl{z) 



q 











P' 





-P' 





r' 





-r' 






-{gKi{z) + Kc(z))l, 



(18) 



where 



76 — (3c 
7C-I-/35 



(19) 
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For the singular case g = 0, we adopte a — 0, p' — j and 
/ = 13. 

The meaning of the new Stokes reference system is as 
follows: after the a rotation, the new generalized Stokes 
parameters Q2 and V2, projections of Ji on axis Q2 and 
V2 correspond to elliptic polarizations, Zeeman elliptic 
absorption affects Q2 only (absorption along axis Q2)- 
Faraday rotation still affects V2 (and U2) but with zero 
absorption. Parameters I2 and U2 are unchanged (the cor- 
respoding axis I and U are not affected by the a rotation) . 
Note that a is wavelength dependent and therefore the ro- 
tation in the Minkowski space is not constant with A! 

In the real space, the meaning of this rotation is not 
any longer as simple as before. However, note that the 
most general state of polarisation is elliptic!! At each A we 
can determine the ellipse of polarisation absorbed by the 
Zeeman effect. We then choose it as axis Q2. The com- 
plete new generalized Stokes system follows from trans- 
formation matrix R2. In deriving a, q — \/\p- + stands 
for the total intensity of the elipse of polarisation and 
sin a = ^f^2^^2 ^^e rate of circular polarisation. Al- 
though easy to calculate, a has no simple meaning in 
terms of the magnetic field, except for the case of a 
strong field when all the Zeeman components are com- 
pletely separated. Then a = for the tt component, and 
tana = zfc l^^fg for the a components, where 9 is the in- 
clination angle of the magnetic field. 

We repeat here all the steps made for first rotation Ri, 
defining the new transformed Stokes and emission vectors 



I2 — R2I1 
J2 = R2J1, 

and calculating the term 



R2) R2"' 



0^ 
da\ [0001 
dz) \ 

0-10 0, 



(20) 



(21) 



so that we can write the transformed transfer equation as 



—I2 = -K2/2 + J2, 
dz 



(22) 



where 



with 



K2 = Ki{z)N2 + (gniiz) + Kciz))t, 













P' 











s 


-r' 


/ 



where we have introduced a new parameter 

da 1 

dz Kl 

This parameter s, a new non-zero entry in K2, makes it a 
little more complicate than before. A new transformation 



is necessary if we want to obtain a simpler matrix like K'2 . 
The way for this simplification is a third rotation R3 given 

by 



R3 = 

with 
cos^ 



10 
1 



COS ^ — sin ^ 
>0 sin^ cos^ 



sin^ 



(23) 



(24) 



The same mathematical steps of the two previous rota- 
tions are repeated for R3. We pass directly to the final 
expression for the transfer equation: 



— /3 =-K3/3 + J3 

dz 

where K3 has the following aspect: 
K3 = Ki{z)N3 + (gmiz) 



(25) 



ith 



^3 = 



/O q 

q 

-p 

Vo 





p 





Kc{z))l, 





r 

0/ 



where we have defined 



p'^ + 
d^l_ 

dz Kl 



(26) 



By now the matrix, and consequently the RTE, has 
been simplified in a general way, without any assumption 
nor constraint. To proceed to an exponential solution for 
the transfer equation, we need to ensure that commuta- 
tion condition (||) holds for K3. A necessary and sufficient 
condition for that is a matrix N3 of the form 



N3 



(N3)o-/(^) 



(27) 



where (N3)o is a constant matrix, and f{z) is any scalar 
function of z. Let P,Q and R be the integrals of p,q and r. 
When calculating the commutator of N3 with its integral 
[N3, L3], its only a priori non-zero entries are (Pq — Qp) 
or {Pr — Rp) . It is very easy to see that these expressions 
vanish when 1) p = 0, 2) g = r = or 3) p^q and r are 
all proportional to the same scalar function f{z). While 
the treatement of cases 1) and 2) is straightforward, the 
general case 3) needs some discussion; we rewrite p,q r 
as pof{z), qofiz) and rof{z). Matrix (N3)o keeps the ap- 
pearance of N3 but with p, q and r substituted by po:9o 
and tq. These new variables to be constant is the nec- 
essary and sufficient condition for writing an expo- 
nential scalar-like solution. A constant K matrix implies 
7 constant parameters. Matrix N3 contains only three pa- 
rameters q,p and r, but we need to keep constant only 



6 



M. Semel & A. Lopez Ariste: The exponential solution 



the two ratios p/q and r/q to satisfy condition (5). Af- 
ter transformation, only two variables are requested to be 
constant, i.e., two constraints instead of seven originally: 5 
degrees of freedom have been earned. Reviewing the trans- 
formation process, those 5 degrees of freedom may be used 
to treat analytically gradients with depth in azimuth, Kc 
and K/, angles a and ^, and f{z). Details about how to do 
it and its application to a numerical code, are left for a 
forthcoming paper. 

For the sake of demonstration, we discuss the variation 
of angles a and 9 alone in the case of separated Zeeman 
components. In absence of azimuth variation p' ~ Q and 
p = s. 

For the tt Zeeman component, polarisation is purely 
linear, a = 0, and we adopt £^ — Q p — s ^ Q: the integra- 
tion is straightforward (case 1). For the a Zeeman com- 
ponents, polarisation is elliptic. Since p' — 0, we adopt 
^ = 7r/4 and r = r' . 

da d(cos 9) 2 

dz A dz 1 + cos^ 9 

Absorption in each a component is proportional to q cx H- 
cos^ 9 and also r = r' cx 1+cos^ 9, 9 being the only variable 
depending on depth. We now suggest f{z) = 1 + cos^ 9 to 
match the depth variation of q, r and P — s — ^-^^ (case 

3) with ^ = 7r/4 To keep the same depth variation for 
we impose a variation of 9 such that 



which is a biquadratic equation, so first we solve 



d{cos9) 



1 



dz {1 + cos^ 



constant, 



and we obtain 
cos 6* 



1 -I- COS"* 



arctan(cos0) = 2 constant(z — zq) 



In general, it will be impossible to interpret a in such a 
simple way. 

3. Diagonalization of the off-diagonal matrix 

In the last section, matrix K was simplified as much as 
possible. All the way to the final form we have seen how 
to incorporate some atmospheric gradients into the ex- 
ponential solution, and at the end we have required the 
commutation condition to still hold with minimum free- 
dom restrictions. We can then calculate the matrix in the 
exponent of the solution; but we must compute the expo- 
nential of this matrix too. The last can be done by using 
the matrix series (^, but we prefer the comfort of cal- 
culating the exponential of scalars. We can achieve this 
purpose by diagonalizing K3. In fact, we only need to di- 
agonalize (N3)o. We solve the eigenvalue equations for this 
matrix: 

{Ns)oU, = X,U, 
The four eigenvalues Xi are given by: 



'0 



where 

and the four eigenvalues are therefore: 



Ai = -V(9o-^o-P§ + A)/2, 

•^2 = — Ai, 

A3 = -^{q^,-rl-pl-A)/2, 
A4 = —A3. 

We introduce the notations: 

9i — — A7; 



9o> 



(29) 



(30) 



(31) 
(32) 

(33) 
(34) 



(35) 
(36) 



noting that 

9i = —92 
83 — —6*4, 

and that 

h = h 



Using this notation, one finds that the corresponding 
eigenvectors Ui are given by 



/ 



\/2go 



go 



Si/pa 
\roSi/{pa9i) , 



(37) 



We may now calculate matrices T and T ■'■ which diago- 
nalize (N3)o: 



and 



T = 




9oPo goPo 9oPo 

P<dQ\ -P0O3 P0O3 

61 53 S3 

-ro6i/9i roS3/93 -roS3/93, 



•-(53 9183/ qo qoPo 

-S3 -6*1^3/50 qoPo 

A\/2 [ Si -93Si/qo -goPo 

<5i 93Si/qQ -qoPo 



hdlpo/{roqo) 



~9i9lpo/{roqo) , , 
'030lpo/{roqo) ' 
939lpo/{roqo) 



Xt + XUrl+pl~q^o)-Qyo=0 



(28) 



These transformation matrices T and T"-"- just derived can 
be applied in general, save for those exceptions where they 
become singular. These particular cases are: 
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1. The case when pq — 0. 

To solve this singularity we suggest the following trans- 
formation matrix and its inverse as a substitution of 
the previous ones: 

-1 0^ 

1 

1 i 

i I, 



and its inverse, so that we get 



/I 
1 



Vo 



(40) 



with the resultant diagonalization of (N3)o 



T(N3)oT- 

















+qo 














-iro 

















(41) 



The case when qo — 0. 

With no loss of generality one can take a = and then 
R2 (and subsequently R3) becomes the identity matrix. 
For the diagonalization transformation, the eigenval- 
ues are the following: 

(42) 
(43) 
(44) 
(45) 



(46) 




(47) 



3. The case when ro = 0, and Qq ^ Pq- 

In this case, A = — Pq ; A2 is imaginary if Qq < Pq. 



The eigenvalues are given as follows: 



Ai 
A2 
A3 
A4 



= iA2, 
= -iA5, 
= 0, 

= 0, 

and the transformation matrix 

/ an «A^ Pq 

Po 
'2qo 







iA2 













2A / 



(48) 
(49) 
(50) 
(51) 



(52) 



The case when ro = 0, and po — qo- 

This is the only case where matrix (N3)o cannot be 

diagonalized. It will be reduced by using 



(53) 




T(N3)oT-^ = 



'0 go 0^ 

go 



>0 0. 



(54) 



Note that the resulting matrix, although not diagonal 
is simple enough as to have its exponential calculated 
as 



J TN3T" 




(55) 



where we have used the notation Q = J qdt. 

Back to the general case, we diagonalize (or apply the 
equivalent transformation in case 4) to the RTE . It is a 
new transformation of the equation, as were Ri, R2 and 
R3, and as in these 3 rotations, a term of the form 



az 



appears in the transformed equation. However, as all the 
entries of (N3)o, and subsequently of T, are constant, this 
term is immediately zero, in perfect agreement with the 
statement that expression (^7|) is a necessary and sufficient 
condition to write a scalar-like exponential solution^ 

For the general case, we define the new generalized 
Stokes and emission vectors after the diagonalizing trans- 
formation for the RTE as 



-Mt + Jt 



(56) 



(57) 



It - TIa 
Jt = T J3 

and rewrite the diagonalized transfer equation as 

dz 
where 

A = Ki{z)T{N3)QT-'f{z) + {gni{z) + k,{z))^ 

is the diagonalized matrix. For the last singular case, the 
equivalent reduced matrix should be used. 

4. The formal solution in a particular atmosphere 

In the previous two sections, we have seen how to trans- 
form matrix K in order to be able to write an exponential 
solution in a convenient form, and to be able to calcu- 
late that exponential as easily as the scalar case. We have 



^ If we relax condition (^^, N3 can still be diagonalized by 
substituting p, q and r for po,qo and ro in T and T^^. But 
this time is not zero. We cannot get rid of this term by 
mathematical manipulations, the RTE cannot be diagonalized 
and a scalar-like solution as Eq. (Kl) cannot be applied. 
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shown that this procedure applies to absorption matrices 
whose depth dependence has a few degrees of freedom but 
still satisfying condition (5). Summing up the last two sec- 
tions, we can manage gradients in the azimuth cj) of the 
magnetic field and in the angles a and ^, variations in k; 
and Kc, treat analytically a general function /(z) in the 
final matrix N3, and integrate for any emission vector J . 

In this context we assume an ideal atmosphere satisfy- 
ing condition (|2^). We write a formal solution (^) for the 
diagonalized transfer equation (p7|): 



•It JT(t)dt + e •''0 It[zo) 



(58) 



Except for the last singular case, we can apply a diagonal 
matrix A, and write its elements as 



KlXif{z) + {gKi + Kc) 



(59) 



where the A^'s coincide with the eigenvalues of the off- 
diagonal matrix (N3)o. We can now write a scalar disen- 
tangled solution for each one of the four components of 
the generalized Stokes vector as 

Jzq 

+e {It)i{zo) (60) 

where we can define four generalized optical depths 

Ti(z,Zo)= / {Kl\if{t') + gKi + Kc)dt', 

J Zo 

complex in general, with which the solutions are 
(ItUz)^ r e^^^'''\jTUt)dt + e^^^'''"\lTUzo). (61) 



We can relate now the generalized Stokes vector It with 
the physical Stokes vector Iq by writting in order all the 
transformations that have been made, i.e. 



It — TR3R2R1/0 — TjIq 



(62) 



where Tj is the complete transformation, product of all 
others in the proper order. Evidently Tj is invertible, so 
that we can make back way from the generalized Stokes 
vector to the physical one: 

lo = Tj'It- 

5. An approximation for a general atmosphere 

Model atmospheres are usually not so perfect as to be 
considered in the cases treated in the last section. In the 
absence of a general analytical solution for the evolution 
operator, we must manage those general models numeri- 
cally. Our strategy is to integrate everything we can and 
to linearize the rest. For this purpose we borrow a tech- 
nique from other successful numerical integrators as the 



well-known DELO (Rees et al. 1989). In what follows we 
develop this idea. 

We want to integrate the transfer equation 



— J= K/ + J 

dz 

or alternatively 

Tz'—^^'^'^- 



(63) 



(64) 



The integration is to be made in the interval zi < z < Z2 
and the atmosphere in the two extreme points zi and Z2 
is specified, that is, we know K(zi) and K(22) and also 
S{zi) and S{z2). The incoming light I{zi) is also given. 
We want to obtain the polarized light at Z2- I{z2)- Af- 
ter the last section we already know how to integrate the 
transfer equation if the atmosphere is characterized by the 
prescriptions given there. Note that, given the atmosphere 
at the two points zi and Z2, one can always calculate the 
angles a and ^ and the matrix N3 at those two levels, and 
look for a suitable integrable atmosphere in agreement 
with expression (27) which satisfy the data at zi and Z2 
as far as possible. So let us approximate our atmosphere 
between these two levels by this integrable atmosphere, 
represented by a matrix K and an emission vector S. We 
can obtain a solution /(Z2) for the equation 



dz 



-K{I-S) 



(65) 



taking as initial condition I{zi) — I{zi). Substraction of 
Eq. (|65|) from Eq. (|3) results in 



-(/-7) = K(7-S)-K(7-S), 

dz 

and upon formal integration : 

I2 72= f \k{I - S) -K{I -S)) dt, 



(66) 



(67) 



where I2 = I{z2) and 72 = I{z2) given as solution to 
Eq. (|65|). This equation reflects the error made under the 
previous approximation. Following our strategy, once we 
have solved for the integrable part we linearize the rest. So 
that we now assume that the right hand side of equation 
( |66| ) is small and can be linearized in the interval zi < z < 
Z2. We define 



y = K{I-S)-K{I-S) 

and upon linearization we write 

y = a + b{z - zi), 
where 



(68) 



(69) 



a = Ki(7i-S'i)-Ki(/i-5i). (70) 
At z — Z2, 'we have y = a + b{z2 — ^i), and we obtain 
a + h{z2 - zi) - K2(72 - ^2) - K2(72 - S2). (71) 



M. Semel & A. Lopez Ariste: The exponential solution 



9 



To solve Eq. m% we write 



1 2 - -f 2 = / ydt. 

'21 



(72) 



And by means of the linearization the last integral be- 
comes 

a{z2 - Zi) + -{Z2 - Zl)^ 

so that, substituting a and b by its complete expressions 

T Z2- Zl 



(K2(/2 - S2) - K2(/2 - ^2 

(Ki-Ki)/i-Ki5i + Ki^i). 



(73) 



In this expression everything is already known except for 
I2 that is precisely what we want to calculate. 

A convenient choice of the overlined parameters may 
render equation ([73[) simpler. For illustration, let us choose 



Ki = K2 = K2 
and Si = Si, S2 ^ 82- We then obtain 



4. The solution being analytical, it is valid for quite thick 
optical layers. The real numerical application is beyond 
the scope of this paper. 

When the commutation condition does not hold, one 
can turn to Magnus' solution, described shortly in the ap- 
pendix. Direct application of Magnus' solution to a nu- 
merical code seems immature at present. For a general 
atmosphere, the numerical strategy proposed is to inte- 
grate analytically what we can and approximate the rest, 
that is: 

1. Divide the atmosphere into a reasonably number of 
layers, so that in each of them the commutation con- 
dition is only slightly violated. 

2. Approximate the general absorption matrix in each 
layer by an average that satisfies the commutation con- 
dition. 

3. Apply the solution developed in this paper using last 
matrix. 

4. Applying an approximation for t he residual matr ix, 
eventually the one used in DELO (Rees et al. 1989). 



12 = 



1 - 



Z2 - Zi 



K2 



Z2 - Zi 



(K2 - Ki) (Ji-Si),(74) 



a solution for I{z2). This solution is not exact, its precision 
depends on how good the linear approximation is. In the 
limit, we can made the integration interval {zi , Z2) as small 
as we want but at the cost of increasing the number of 
layers. A compromise will be necessary between speed and 
required precision. 

6. Conclusions 

The purpose of this paper was first to deepen our un- 
derstanding of the integration of the RTE for polarised 
light and next to improve the basis for numerical codes. 
The main conclusions is: The fundamental key to solve the 
RTE of polarised light is the commutation of the absorp- 
tion matrix and its integral. 

When this commutation condition is satisfied : 

1. A scalar like solution can be proposed to the vector 
equation. 

2. A constant absorption matrix satisfies the commuta- 
tion requirement, however, it is only a sufficient con- 
dition, not necessary. After some elaboration one can 
show that only two constraints at all are necessary in- 
stead of the seven inherent in the fully constant ab- 
sorption matrix. 

3. In general, it will possible to diagonalize the absorp- 
tion matrix and consequently also the RTE with its 
vector solution. This results in four scalar equations 
with four scalar solutions. The variables are no more 
the usual Stokes parameters, but generalized ones, cor- 
responding to general states of polarisation. 



As an objective, we intend to improve the efficiency 
bf integration and inversion codes. This will be a must 
in treating the abundant data expected from multi-line 
spectropolarimetric observations to be provided by the 
French— Itahan telescope THEMIS. 

Appendix A: The general solution to the 
evolution operator 

We intend in this appendix to give a self-contained proof 
of Magnus' exponential solution to the equation for the 
evolution operator. The original proof is to be found in 
Magnus' paper(1954). Here we will give a short but com- 
plete proof introducing in the meanwhile the necessary 
tools used while working with Magnus' expressions and to 
get acquainted with algebraic manipulations inherent to 
problems involving non-commutative operators. Magnus 
makes use of an original technique to manage the deriva- 
tive of the exponential of a linear operator (a square ma- 
trix in our case). He transforms this derivative into an 
algebraical expression, and uses it to solve the differential 
equation. 

In what follows w, x, y and z stand for such opera- 
tors also called Lie elements. The Lie product of two Lie 
elements x and y is defined by 

w = [x, y] , 

where w is a new Lie element. This product is usually 
called commutator. Now an abreviated notation is used 
for the 1-fold Lie product by a; of y as 



times 



{y,A = [[■■■[y,- 
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with {y, x°} = y. With this notation, it is easy to show by 
a straightforward calculation that 



OO 



(A.l) 



1=0 



where we remind that the exponential of an operator x is 
defined as 

OO _j 

Following Magnus, we extend the previous notation to a 
polynom P{x) in an evident form: 



{y,P(x)} = ^Pi{y,x'} 
where 



(A.2) 



1=0 



1=0 



And we are ready to demonstrate the first two formulas 
which we will use hereafter. 
Formula 1 



d 



e^ - 1 



This formula arises from the straightforward calculation 
of its left-hand-side. We begin calculating the effect of y-^ 
on the left of the exponential: 

(^^) '^'^ '"V'^\^y^^^y)^\\ (yx"^ + xyx + x'^y) + . . . 
Next we multiply on the left by e~^: 

= y + ^{yx - xy) + -^^ {x^y - 2xyx + yx^) + . . . = 



11 ^ 

= y + 77 [y, a;] + ^ [[y, + ■■■ = Y1 TTTTv''-^' 



3! 



Now, 



OO -j^ 

^ (T+T)! 

1=0 ^ ' 



e^ - 1 



so Formula 1 is demonstrated. 
Formula 2 



d \ \ _ f 1 - e"^ 



The left-hand-side of Formula 2 is obtained by multiplying 
the left-hand-side of Formula 1 on the left by e^ and on 
the right by e~^ . If we do the same multiplications in the 
right-hand-side of Formula 1 we obtain 



or expanding the right-hand-side 



OO OO 



1=0 n=l ^ ' ^ ' n=0 



where 



(-1) 



which can be seen to correspond to the expansion of 

1 - e-^ ■ 

y, 

X 

Next, Magnus demonstrates what can be called 

The Magnus' Inversion Lemma: Let P{x) and 
Q{x) be two power series in x which satisfy 

P{x)Q{x) = 1. 

Then each of the equations 

{y, Pix)} = y ^ {u, Q{x)} 

is a consequence of the other one. Let P{x) — '^iPix^ and 
Q{^) = J2ni Irax"^, then by hypothesis 



1 = ^Piqn 



Now we can obtain the following equivalent expressions 



{yA}={ y,Y.P^qn.x'^"' \ = 



where we have used notation (A. 2). Next we can separate 
indexes and write 



E'?™|e^'{^'-'}'-'"|-e 



q^{{y,Pix)},x^}. 



If we suppose now that {y, P{x)} = u, immediately we 
obtain that 



The inverse implication is completely equivalent. 

With the Inversion Lemma and Formula 2, we have 
all the instruments to solve the equation for the evolution 
operator 

Exponential solution Theorem (Magnus): Let 

K(i) be a known function of t in an associative ring (for 
our purposes it is a matrix), and let 0{t) be an unknown 
function (in our case the evolution operator) satisfying 



f=KO, 0(0) = 1, 



(A.3) 
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where 1 is the identity matrix. Then, if certain unspecified This expression can be finally expanded using the foUow- 
conditions of convergence are satisfied, 0{t) can he written ing power series 
in the form 

where 



0{t) = expn{t) 



n 



1 



n 



1 - e- 



n=0 



The Pn vanish for n = 3,5,7,..., and P2m = 
(-l)™-iB2m/(2m)!, where the B2m (form = 1,2,3,...; 
are the Bernoulli numbers. 

Integration of this equation by iteration leads to an in- 
finite series for f2 the first terms of which (up to terms 
involving K 'sj are 



where the /3„ has the given values. 

To integrate the resulting equation for Q we start at 
t=0, where ri(O) = to satisfy boundary conditions. In- 
troducing this solution into the equation we obtain a new 
solution: 



KdT. 



We can iterate the procedure to obtain ^2, ^3, 
dflr, 



. as 



dt 



= ^/3„{K,17;;_J. 



n=0 



n{t) 



K(r)dr 



1 



1 
1 
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K(r), 
K(t), / K{a)dcr 



K(t), / K{a)da 
Jo 

Kia), rKip)dp 
Jo 



dr^ 



da 



K{a)da 



1 

24 



K(r), 



K(r) 



K(a), 



K(a) 



Kip), / Ki,.)diy 



dT+ 
dT+ 
,dp 



K{p)dp] }da 



dT- 



Let us suppose 0{t) — expQ{t), then 



dO 

'dt 



dn d 
'dtm 



c-"e". 



By using Formula 2 with y = ^ and x = one obtains 



dO 

~dt 



dn i-c-^^ 



which compared with Eq.(A.3) gives 
^ , , {dn I- c-^^ 



We can now apply the Inversion Lemma with the same 
substitutions in x and y as before, and with u = K and 



1 - e 



to obtain 



and 



n 



dfl 
'dt 



K, 



1 - e- 



The solution is obtained as the limit of this series: 
Q(t) — lim flm- 

m — >oo 
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